w=1;;c=zeros(1,13)';d=[]; S1=[];%b=zeros(1,12)'
while w<10
    m=rand(1,1000000);%Generate a uniformly distributed random sequence between 0 and 1
    M=2.5*(1./(1-m)).^(1/40);%Generate a sequence that follows the Pareto distribution
    M=M-2.5;
    total_sum=sum(M); %The sum of sequence M
    r=sort(M);%Sort the sequence M
    Mmax=max(M);Msum=sum(M);Dmax=(Mmax/2.5)^(1/3);
    %%Find the values in each diameter grade
    v1=find(r>0 & r<=0.001);v2=find(r>0.001 & r<=0.005);v3=find(r>0.005 & r<=0.01);
    v4=find(r>0.01 & r<=0.05);v5=find(r>0.05 & r<=0.075);v6=find(r>0.075 & r<=0.25);
    v7=find(r>0.25 & r<=0.5);v8=find(r>0.5 & r<=1);
    v9=find(r>1 & r<=2);v10=find(r>2 & r<=5);v11=find(r>5 & r<=10);v12=find(r>10 & r<20);
    k=abs(0+1*randn(1,12));%Generate positive numbers following a normal distribution
    %Group the values in each interval
    s1=r(1,v1)*k(1);s2=r(1,v2)*k(2);s3=r(1,v3)*k(3);s4=r(1,v4)*k(4);s5=r(1,v5)*k(5);
    s6=r(1,v6)*k(6);s7=r(1,v7)*k(7);s8=r(1,v8)*k(8);s9=r(1,v9)*k(9);s10=r(1,v10)*k(10);
    s11=r(1,v11)*k(11);s12=r(1,v12)*k(12);
    %s1=s1.^3;s2=s2.^3;s3=s3.^3;s4=s4.^3;s5=s5.^3;s6=s6.^3;s7=s7.^3;s8=s8.^3;s9=s9.^3;s10=s10.^3;s11=s11.^3;
    %Calculating the soil weight percentage in each grade
    s=[sum(s1),sum(s2),sum(s3),sum(s4),sum(s5),sum(s6),sum(s7),sum(s8),sum(s9),sum(s10),sum(s11),sum(s12)];
    q=sum(s);
    B=[s(1)/sum(q),s(2)/sum(q),s(3)/sum(q),s(4)/sum(q),...
        s(5)/sum(q),s(6)/sum(q),s(7)/sum(q),s(8)/sum(q),...
        s(9)/sum(q),s(10)/sum(q),s(11)/sum(q),...
        s(12)/sum(q)]'; b1=flipud(B);
    Y=cumsum(B).*100; 
    X=[20;10;5;2;1;0.5;0.25;0.075;0.05;0.01;0.005;0.001];
    f=inline('a(1)*x.^(-a(2)).*exp(-x/a(3))','a','x');
    [aa,resnorm]=lsqcurvefit(f,[100,0.01,10],X,Y);%calculating the parameters
    y=aa(1)*X.^(-aa(2)).*exp(-X/aa(3));%examine the R-suqared
    SSE=(Y-y).^2;
    SSR=(y-mean(Y)).^2;
    SST=sum(SSE)+sum(SSR);
    RR=sum(SSR)/SST;
    if  RR>0.99 %constraints
        d(:,w)=B;
        A(:,w)=aa;%A(4,w)=RR;
        %draw a frequency distribution graph
        pos1 = [0.1 0.1 0.3 0.3];
        subplot('Position',pos1)
        z=[0.001,0.005,0.01,0.05,0.075,0.25,0.5,1,2,5,10,20];
        z1 =0.001:0.0001:20;
        n1 = interp1(z,b1,z1,'pchip');
        semilogx(z,b1,'*',z1,n1)
        set(gca,'ytick',0:0.1:0.8);
        xlabel('D');ylabel('P(D)(%)');title('Frequency distribution curve');
        hold on
        %draw GSD frequency curve
        pos2 = [0.5 0.15 0.4 0.7];
        subplot('Position',pos2)
        o=fliplr(z)';
        xx=o./aa(3);%calculating D/Dc
        yy1=100*((o.^aa(2)).*cumsum(B))./aa(1);%the GSD function
        yy=yy1./max(yy1);
        loglog(xx,yy,'*');
        title('GSD fitting curve');
        hold on
        %draw cumulative distribution curve figure
        pos3 = [0.1 0.6 0.3 0.3];
        subplot('Position',pos3);
        c(1)=1;
        c(2:13)=1-cumsum(B);
        S1(:,w)=c;
        z1=0.0001:0.0001:20;
        XX=[20;10;5;2;1;0.5;0.25;0.075;0.05;0.01;0.005;0.001;0.0001];
        c1=interp1(XX,c,z1,'pchip');
        semilogx(XX,c,'*',z1,c1)
        set(gca,'XDir','reverse')
        axis( [0.0001 20 0 1] )
        title('cumulative distribution lines');%add title
        grid on
        hold on
        w=w+1
    else
        w=w
    end
end
subplot('Position',pos2);
x=0.0001:0.1:20;
y=exp(-1.016*x);
loglog(x,y);
ylim([0.01 1]);
hold on





